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We simulate the motion of spherical particles in a phase-separating binary mixture. By combining 
cell dynamical equations with Langevin dynamics for particles, we show that the addition of hard 
particles significantly changes both the speed and the morphology of the phase separation. At 
the late stage of the spinodal decomposition process, particles significantly slow down the domain 
growth, in qualitative agreement with earlier experimental data. 
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. Phase separation plays a significant role in determining the morphology and properties of polymer composites, 
which typically involve a blend of various macromolecular fluids and solid "filler" particles [|J. Despite the utility of 
these composites, there is little understanding of the kinetic processes (including phase separation and wetting) that 
occur in the complex mixtures. While phase separation in binary systems has been studied extensively theoretically 
and experimentally ||||, the influence of solid additives on the mixtures is still poorly understood. Recent studies 
have shed light on the interactions between a phase-separating fluid and a stationary wall Gj , sphere Q or substrate 
but much less is known about the kinetics of mixtures that contain mobile particles. To address this problem, 
Tanaka et al. || examined the properties of a polymeric mixture undergoing a critical quench in the presence of 
small glass particles, which are preferentially wet by one of the components. Their results revealed that even a small 
concentration of hard particles significantly changes the morphology and dynamics of the phase separation process. 
However, no theoretical or computational model was developed to characterize these changes. 

In this Letter, we report the first simulations of hard mobile particles in a phase-separating binary mixture. Unlike 
earlier dynamical models of ternary systems (developed mostly for oil- water-surfactant mixtures [p|-|l2| ) , we explicitly 
take into account the "excluded volume" interaction between the particles and the background fluid. Furthermore, 
we can vary the particle-fluid interactions, allowing for a richer range of behavior than that of a surfactant. Thus, the 
model presents a new means of exploring the physical properties of complex mixtures containing colloidal particles. 
Here, we consider particles that are preferentially wet by one of the two components and show that the boundary 
00 ' and "excluded volume" conditions at the particle surfaces significantly slow down the domain growth and change the 
£NJ \ morphology at the late stage of the phase separation. 

We consider a phase-separating symmetric binary AB mixture that is characterized by the scalar order parameter 
| 'J. The phase separation dynamics are described by the Cahn-Hilliard equation, 

where T is a kinetic coefficient, £ is a conserved zero mean Gaussian white noise with covariance (£(r, f)£(r', t')) = 
— GiV 2 <5(r — r')S(t — t'), and T is a free-energy usually given by the Ginzburg-Landau functional, 



T= |dr{-^ 2 + !* 4 + |(W) 2 }. (2) 

Into this system, we introduce spherical particles of radius Rq that undergo Brownian motion. The particle dynamics 
are described by the following Langevin equation, 
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c3 : R i = Mf i +r ?i) (3) 

where M is the mobility, fi is the force acting on the i-th particle due to all the other particles, and r/ represents 
Gaussian white noise with (r)i a (r, t)r)jp(r', t')) = G2^(r — r')S(t — t r )5ij5 a p. In this study, we neglect interactions 
between particles (i.e., fi = 0) and only take into account the particles' diffusive motion. We also disregard osmotic 
effects (i.e., coupling between the particle motion and the order parameter field). 

The simulation is carried out in two dimensions; our lattice is 256 x 256 sites in size, with periodic boundary 
conditions in both the x and y directions. A cell dynamical systems (CDS) method Jll| is used in place of a direct 
forward integration of Eq. ([!]) to update the value of ^ for the phase-separating AB mixture. Note that ^ — 1 (-1) 
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corresponds to the equilibrium order parameter for the A-rich (B-rich) phase. By employing CDS modeling (rather 
than a conventional discretization of Eq. (|l|)), we can significantly increase the computational speed of the simulation. 
To simulate the particle dynamics, we discretize Eq. (||) and only allow the particles to move between different lattice 
sites. A "Kawasaki exchange" mechanism is used for each particle move: first, the order parameter values from all 
the cells to be occupied by a particle in its "new" position are moved to the "old" position; next, the boundary and 
excluded volume conditions are imposed for the order parameter at the "new" particle position. This mechanism 
ensures the conservation of the order parameter. Such dynamics may break down for high particle mobilities, so we 
considered only the case where the diffusion constant is rather low (almost all particle "jumps" are to neighboring 
sites). The discretized equations of motion have the following form, 

*(r, t + 1) = F[*(r, t)] - {{F[9(t, t)] - *(r, i)» + £(r, t), 
F[¥(r, t)} = /(¥(r, t)) + D(«*(r, t)» - *(r, t)), 

/(*) = Atanh(tf), 
R s (t+l)=R i (t) + Mf i + r ?i (t), (4) 

where ((*)) is the isotropic spatial average over the nearest-neighbor and the next-nearest neighbor sites, and [((*)) — *] 
can be thought of as a discrete generalization of the Laplacian. 

At the surface of each particle, the lattice boundary conditions (specified order-parameter value and zero order- 
parameter flux) are imposed as: \&(r, i) = fy s , and d n F(r,t) = 0, if Rq < |r — Rj(i)| < Ro + a , where a is the lattice 
spacing and d n denotes the "lattice" normal derivative. Here, we set $ s = 1 so that the particles are "coated" by 
fluid A. The d n F = condition ensures zero flux of \P into the particles since F plays the role of a chemical potential. 

The function F in Eq. (^J) has a local driving term / and a term arising from the interaction with other sites; the 
map controls the local dynamics at each site. It is critical that / has a single unstable fixed point and two stable 
fixed points symmetrically located on each side of the unstable fixed point. Its exact functional form is not important 
for studing the universal properties of the phase separation dynamics jl3| . Here, we select the map /(^) = Atanh(W), 
with A < 1 above the critical temperature, and A > 1 below. 

We perform simulations for systems containing 0, 25, 50, 100, 150, 300, and 400 particles of radius 1 (all lengths 
are given in units of the lattice spacing a) |Q. Each system was averaged over 3 runs of 20,000 time steps each. 
For all systems, the following values of the parameters were used: A — 1.3, D = 0.5, G\ = 0, Gi = 0.5 . The initial 
fluctuations of \& are Gaussian with a variance of 0.05. For all runs, the composition of the fluid is fixed at 50:50 
(representing a critical mixture). 

The characteristic length R(t) of the evolving domains is plotted as a function of time in Fig. |I|a. For R(t) we use 
the "broken bond" formula [fl5j| , R ~ L d /A(t), where L d is the volume of the system and Ait) is the total interfacial 
"area" . For d — 2 this becomes, 

T 2 

R= , (5) 

N X + Ny ' ^ > 

where L is the system size, N x and N y are the numbers of "broken bonds" (pairs of nearest neighbors with opposite 
signs of \& in the x and y-directions, respectively). This measure of a characteristic length empirically yields the 
correct asymptotic behavior for both critical and off-critical quenches in binary mixtures. 

The simulations reveal that the presence of particles slows down the domain growth in the late stage. It can be 
seen that for large particle numbers (N > 100), R(t) undergoes a change from a Lifshitz-Slyozov jl6| regime, with the 
growth exponent of 1/3, to a new regime. This new, slow-growth behavior is also characterized by a new morphology, 
different from the critical or slightly-off-critical pattern of bicontinuous domains. p7[ The change in morphology can 
be seen by comparing order parameter patterns and particle positions for the system N — 300 at times t=300 ( Fig. 
^|a) and t=3000 (Fig. |^b). The slower domain growth and altered structure of the mixture are qualitatively similar 
to the observations of Tanaka et al. || 

In the early stage of phase separation, the formation of interfaces occurs as in a normal critical mixture, as would 
be expected as long as the concentration of particles is small, i.e., the interparticle distance is much larger than the 
particle radius. During this time, the initial domain growth satisfies the Lifshitz-Slyozov law, with the prefactor E(n) 
smoothly depending on the particle density, 

R{t) = E^t 1 / 3 , (6) 

where n = N/L 2 is the particle density, and E(n) = Eo(l + an), with Eq w 0.4 being a growth prefactor in a 
particle-free system. This growth continues until the characteristic domain size becomes comparable to the average 
interparticle distance rT 1 ! 2 . At that point, the wetting phase (A) percolates to form a single infinite domain (see Fig. 
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, and droplets of the B-phase are trapped inside this domain. Further coarsening of B-domains is inhibited by the 
particles acting as obstacles to the motion of interfaces. 

To describe the dependence of the characteristic size on both time and the particle density, we use the following 
scaling function, 

R(t) = rT l ' 2 {l + an)G{btn^ 2 ), (7) 

where b « 1 is a "metrical factor", G(x) behaves as G{x) ss Gqx 1 ^ 3 for small x (6 1 / 3 Go = Eq), and 7 = 3 (which 
is required to satisfy the transition to the Lifshitz-Slyozov growth law for n — > 0). There is clearly a slowing of 
growth at large time, and it is reasonable to assume G(x) ~ x b with a small power S (or even logarithmic growth) 
for large x. A similar scaling form as eq. (^) was used by Gyure et al. [Gl| to describe the dependence of the domain 
growth on the number of impurities in an Ising model. We introduce the additional factor 1 + an to account for the 
effective off-criticality induced by particles. To illustrate this scaling behavior, we plot the characteristic size in scaled 
coordinates p = R(t)n 1 ^ 2 / (1 + an) vs. r = tn 3 / 2 . For a = 13.1, it can be seen (Fig. |l]b) that all data fit reasonably 
well onto one master curve, with the exception of the N — 400 case, where additional n-dependence is presumably 
required. 

The observed slowing down of the domain growth is reminiscent of the interface pinning in Ising-type systems 
with quenched impurities p8|-^l|. In all those studies, impurities reduced local interfacial tension and thus enforced 
late-time pinning, with domain growth slowing down logarithmically, R oc (lnt) 77 . A similar effect was also seen in the 
"hybrid" model of Kawakatsu et al. [jl2j for a binary mixture with surfactant molecules. On the other hand, our hard 
particles, which do not behave as surfactants and prefer to be in the bulk (A) phase, act as obstacles to the interface 
motion. When the characteristic domain size becomes comparable to the interparticlc distance, interface coarsening 
becomes hindered by these obstacles, and the slowing down occurs. It is likely that at the very late stage, the domain 
growth would stop completely (as indicated by the behavior of the N = 400 curve); however, it is difficult to verify 
this hypothesis computationally. 

To elucidate the importance of wetting on the slowing down, we performed a simulation with hard, mobile particles 
and no preferential adsorption, i.e., with the following boundary conditions on the surface of the particles: <9„'I'(r, t) — 
0,d n F(r,t) = 0. For N=300 particles and no wetting, we observed no deviation from the Lifshitz-Slyozov growth law 
within the time scale of our simulations (20,000 timesteps) . This result indicates that the slowing down is clearly 
enhanced by the strong wetting and not merely by the effects of excluded volume or particle mobility. Indeed, non- 
wetting particles neither pin nor block interfaces, and thus, have only minimal effect on the dynamics of the late-stage 
coarsening. 

Although the simulations were performed for two-dimensional systems, we believe that the major features (the 
initial off-criticality and the late-stage slowing down) would be found in three dimensions as well. Indeed, in 3d, 
the particles still would represent a network of obstacles for the coarsening 2d-interfaces. However, verifying this 
prediction remains a major computational challenge, even though the proposed model can be easily extended to three 
dimensions. 

In the computations described here, we have not taken into account hydrodynamic interactions or the dependence of 
viscosity on the order parameter. We also have not considered systematically the dependence of the growth behavior 
on such factors as temperature (which manifests itself in the effective diffusion constant D e ff), the clustering of 
particles that could lead to additional domain pinning (an effect observed by Tanaka et al. || at very high particle 
densities), the dependence on the nature of the interparticlc potential and other features. These phenomena will be 
the subjects of future studies. 
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FIG. 1. Characteristic length as a function of time: a) unsealed R vs. t, b) scaled coordinates p = Rn 1 ^ 2 /(l + an) vs. 
r = tn 3 ^ 2 , where n = N/L 2 (data averaged over 3 runs). Growth exponent v ~ 1/3 in the early time range 1CT 2 < r < 10°. 

FIG. 2. Order parameter patterns and particle positions for N=300 particles at t=300 (a) and t=3000(b). Black circles 
represent particles, dark grey regions are the A-phase domains, light grey regions are the B-phase domains, and white points 
are interfaces ft* 0). 
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